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A finite size scaling is applied to the Yang-Lee zeros of the grand canonical partition function 
for the 2-D Hubbard model in the complex chemical potential plane. The logarithmic scaling of 
the imaginary part of the zeros with the system size indicates a singular dependence of the carrier 
density on the chemical potential. Our analysis points to a second-order phase transition with critical 
exponent | = | ± | which leads to a divergence of the electronic susceptibility. We interprete these 
results as reflecting singular behaviour of the density of states in the quasiparticle spectrum. 
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I. INTRODUCTION 

Quantum Monte Carlo techniques are widely applied 
in the investigation of the phase structure of the Hub- 
bard model in an attempt to explore its relationship to 
high temperature superconductivity (HTSC). The pio- 
neering works of Blankenbecler, Scalapino and Sugar Jl]] 
and Hirsch outlined the finite-temperature Quantum 
Monte Carlo algorithms for treating systems of interact- 
ing fermions, thus allowing efficient numerical simula- 
tions of the model. These simulations showed remarkable 
stability for the half-filled model with on-site Coulomb re- 
pulsion, and could be extended to the doped case in the 
relatively high temperature regime. However, in the low- 
temperature and finite-doping regime, which is more rele- 
vant to HTSC, the simulations become less reliable due to 
the frequent appearance of field configurations with neg- 
ative statistical weight, the so-called sign problem. This 
problem severely obscures numerical evaluation of ther- 
modynamical averages of correlation functions and local 
condensates with subsequent difficulties in establishing 
the phase structure of the Hubbard model. 

Barbour and Klepfish |3| explored an alternative 
method of detecting a phase transition in the Hubbard 
model. This method is based on the application of the 
general ideas of the Yang-Lee theorem Q regarding the 
distribution of zeros of the grand canonical partition 
function in the complex fugacity plane. The theorem 
was originally proved for the fsing model and lattice 
gas[Q. Later extensions of this theorem include the case 
of higher-order Ising spins Ising models with multiple 
spin interactions, the quantum Heisenberg model the 
classical XY and Heisenberg models and some con- 
tinuous spin systems §. Recently, a generalization of the 
Yang-Lee theorem to the asymmetric first-order transi- 
tion was madeQ. 

Following the Yang and Lee treatment, one identifies a 



phase transition by the appearance of zeros of the parti- 
tion function in the complex chemical potential (/i) plane 
near the real axis, namely the physical region. For any 
finite system simulated on a lattice, these zeros will never 
be real. However, one may identify a group of zeros in 
the vicinity of the real axis which, if there is a phase tran- 
sition, converge towards it as the system volume grows. 
The critical value of the chemical potential is thus the 
thermodynamic limit of the zero with the smallest imag- 
inary part. 

The 4 2 lattice simulations presented in Ref.|| sug- 
gested that, from the scaling of the smallest zeros in 
the complex /i plane with the inverse temperature (3, the 
two-dimensional single-band Hubbard model undergoes 
a phase transition only at zero temperature. This con- 
clusion is in principle supported by the Mermin- Wagner 
theorem which states the absence of the off-diagonal long 
range order in two-dimensional systems at temperatures 
above zero (jTo|. However, it is not clear whether this 
theorem excludes any kind of phase transition in a doped 
two-dimensional Hubbard model. Therefore a thorough 
finite size analysis is required to validate or reject the 
above conclusion, as was pointed out in j|. 

The distribution of zeros in the vicinity of the real axis, 
namely the angle of the zeros' locus and the scaling of the 
location of these zeros in the complex plane, allows one 
to derive the critical exponents for the studied system 
flU , Q , |l3| . This derivation is based upon the scaling 
hypothesis which is fundamental to the renormalization- 
group approach. In this work we apply the ideas devel- 
oped in statistical mechanics models to investigate the 
criticalities in the Hubbard model of strongly correlated 
electrons. 

One such criticality manifests itself through the macro- 
scopic density of carriers. When the latter is derived by 
differentiating the free energy with respect to the chem- 
ical potential, it can serve as an indicator of a phase 
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transition controlled by the chemical potential. As in 
order-disorder transitions, one would expect a symmetry 
breaking signalled by an order parameter. In this model, 
the particle-hole symmetry is broken by introducing an 
"external field" which causes the particle density to be- 
come non-zero. Furthermore, the possibility of the free 
energy having a singularity at some finite value of the 
chemical potential is not excluded: in fact it can be a 
transition indicated by a divergence of the correlation 
length. A singularity of the free energy at finite "exter- 
nal field" was found in finite-temperature lattice QCD 
by using the Yang-Lee analysis for the chiral phase tran- 
sition |Q. A possible scenario for such a transition at 
finite chemical potential, is one in which the particle den- 
sity consists of two components derived from the regular 
and singular parts of the free energy. 

Since we are dealing with a grand canonical ensemble, 
the particle number can be calculated for a given chem- 
ical potential as opposed to constraining the chemical 
potential by a fixed particle number. Hence the chem- 
ical potential can be thought of as an external field for 
exploring the behaviour of the free energy. From the mi- 
croscopic point of view, the critical values of the chemical 
potential are associated with singularities of the density 
of states. Transitions related to the singularity of the 
density of states are known as Lifshitz transitions [ p^[ . 
In metals these transitions only take place at zero tem- 
perature, while at finite temperatures the singularities 
are rounded. However, for a small ratio of temperature 
to the deviation from the critical values of the chemical 
potential, the singularity can be traced even at finite tem- 
perature. Lifshitz transitions may result from topological 
changes of the Fermi surface, and may occur inside the 
Brillouin zone as well as on its boundaries fll6fl . In the 
case of strongly correlated electron systems the shape of 
the Fermi surface is indeed affected, which in turn may 
lead to an extension of the Lifshitz-type singularities into 
the finite-temperature regime. 

In relating the macroscopic quantity of the carrier den- 
sity to the density of quasiparticle states, we assumed the 
validity of a single particle excitation picture. Whether 
strong correlations completely distort this description is 
beyond the scope of the current study. However, the iden- 
tification of the criticality using the Yang-Lee analysis, 
remains valid even if collective excitations prevail. 

The paper is organised as follows. In Section 2 we out- 
line the essentials of the computational technique used 
to simulate the grand canonical partition function and 
present its expansion as a polynomial in the fugacity vari- 
able. In Section 3 we present the Yang-Lee zeros of the 
partition function calculated on 6 2 -T0 2 lattices and high- 
light their qualitative differences from the 4 2 lattice. In 
Section 4 we analyse the finite size scaling of the Yang- 
Lee zeros and compare it to the real-space renormaliza- 
tion group prediction for a second-order phase transition. 
Finally, in Section 5 we present a summary of our results 



and an outlook for future work. 



II. SIMULATION ALGORITHM AND FUGACITY 
EXPANSION OF THE GRAND CANONICAL 
PARTITION FUNCTION 

The model we are studying in this work is a two- 
dimensional single-band Hubbard Hamiltonian 



- ^^{n i+ + 7ij_) 



(1) 



where the i, j denote the nearest neighbour spatial lat- 
tice sites, a is the spin degree of freedom and rn a is the 
electron number operator c\ a Ci„. The constants t and 
U correspond to the hopping parameter and the on-site 
Coulomb repulsion respectively. The chemical potential 
/i is introduced such that [i — corresponds to half- 
filling, i.e. the actual chemical potential is shifted from 
fx to (1 — ¥. 

The finite-temperature grand canonical partition func- 
tion is given by 

Z = Tr (e~ Pk ) ■ (2) 

Following Hirsch 0] and White et al. jl8| we rewrite Z 
as 



Tr 



no 
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where K corresponds to the nearest neighbour hopping 
term in the Hubbard hamiltonian (|l|) and V to the on-site 
interaction including quartic products of fermion fields, 
and N to the particle number operator. This decom- 
position, based on the Trotter formula |^(|, introduces 
a systematic error proportional to dt 2 . Thus Z can be 
represented as a path integral of a 2+1-dimensional sys- 
tem where the range of the additional dimension (imag- 
inary time) corresponds to the inverse temperature via 
/3 = dt x n T . The quartic interaction can be re- written 
in terms of Ising fields Si t i using the discrete HS trans- 
formation [|| 



-dtv 
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where i, I is the space-time index of a lattice site and 
the coupling 7 is related to the original on-site repulsion 
constant by 



cosh (dtj) = exp 



(5) 
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This transformation enables one to integrate out the 
fermionic degrees of freedom and the resulting partition 
function is written as an ensemble average of a product 
of two determinants 

Z = * = E det(M+) det(M~~) (6) 

{s M =±l} { Si ,,=±l} 

such that 



m ± = (/+p ± ) = (/ +n s ; d 
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where the matrices Bf 1 are defined as 



B ± = e -(±dtv) e -dtK e dt» 



(7) 



(8) 



with Vij = SijSij and — 1 if i,j are nearest neigh- 
bours and Kij = otherwise. The matrices in (fy and (||) 
are of size (n x n y ) x (n x n y ), corresponding to the spatial 
size of the lattice. 

The expectation value of a physical observable at 
chemical potential fi, < O > /1 , is given by 

„ , . E 0({su})i(p,{s hl }) 

< o > = I M = {s " fc±l} (9) 

{«i,j=±i} 

where the sum over the configurations of Ising fields is 
denoted by an integral. Since z(ji) is not positive definite 
for Re(/i) ^ we weight the ensemble of configurations 
by the absolute value of z(n) at some fx = /io- Thus 



< O > L 
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The partition function Z(p) is given by 



Z(ji) 
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(ii) 



The angular brackets notation (<>>„) is to be under- 
stood here and in what follows as Monte Carlo integra- 
tion with importance sampling based on the absolute 
value of the z calculated at an update chemical potential 
fig. The normalization of Z is irrelevant as can be seen 
from eqn.(|lC|). 

It was shown in |J that a unitary transformation 
equivalent to the particle-hole transformation allows the 
partition function to be represented as a power series in 
e M/3 or e M/3 _|_ e -M/3 w hose powers are ranging between 
[—n x n y , n x n y ] or [0,n x n y ] respectively, where the coeffi- 
cients of the first expansion are the canonical partition 
functions for a number of particles corresponding to the 



index n, with n = being the half-filling partition func- 
tion. These coefficients are obtained by ensemble aver- 
aging of the following expansion of the statistical weight 
for each configuration 



or alternatively 



fm m n y /3 
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where the expansion coefficients in both these equations 
are obtained from the eigenvalues of the matrix P,~ „ 
B . Alternatively, each of these expansions can be done 
around the updating fugacity point, thus yielding for 
eqn.(E|) 



n=0 



(14) 



with fiQ being the updating value of the chemical poten- 
tial. A similar expression can be obtained for eqn.(fl2"|). 
The coefficient of the zeroth power in (|l4|), normalised 
and averaged over the ensemble of the HS field configu- 
rations, is the average sign of the statistical weight cal- 
culated as 



< sign >= 
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where N c is the number of the HS field configurations 
accepted in the Monte Carlo integration. The statisti- 
cal fluctuation in this coefficient reflects the sign fluctu- 
ations. The expansion for the partition function is then 
given by: 



Z(fi) = e 



[in x n y (3 
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When the average sign is near unity, it is safe to as- 
sume that the lattice configurations reflect accurately the 
quantum degrees of freedom. Following Blankenbecler et 
al. ^] the diagonal matrix elements of the equal-time 
Green's operator G ± = (I + P ± )^ 1 accurately describe 
the fermion density on a given configuration. In this 
regime the adiabatic approximation, which is the basis 
of the finite-temperature algorithm, is valid. The situa- 
tion differs strongly when the average sign becomes small. 
We are in this case sampling positive and negative z(/io) 
configurations with almost equal probability since the ac- 
ceptance criterion depends only on the absolute value of 
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In the simulations of the HS fields the situation is dif- 
ferent from the case of fermions interacting with dynam- 
ical boson fields presented in Ref.||]. The auxilary HS 
fields do not have a kinetic energy term in the bosonic 
action which would suppress their rapid fluctuations and 
hence recover the adiabaticity. From the previous sim- 
ulations on a 4 2 lattice || we know that avoiding the 
sign problem, by updating at half-filling, results in high 
uncontrolled fluctuations of the expansion coefficients for 
the statistical weight, thus severely limiting the range of 
validity of the expansion. It is therefore important to 
obtain the partition function for the widest range of (j,q 
and observe the persistence of the hierarchy of the ex- 
pansion coefficients of Z. An error analysis is required 
to establish the Gaussian distribution of the simulated 
observables. We present in the following section results 
of the bootstrap analysis jl7| performed on our data for 
several values of fio- 

III. TEMPERATURE AND LATTICE-SIZE 
DEPENDENCE OF THE YANG-LEE ZEROS 

The simulations were performed in the intermediate 
on-site repulsion regime U — At for f3 = 5, 6, 7.5 on lat- 
tices 4 2 , 6 2 , 8 2 and for (3 — 5, 6 on a 10 2 lattice. The ex- 
pansion coefficients given by eqn. (|l4|) are obtained with 
relatively small errors and exhibit clear Gaussian distri- 
bution over the ensemble. This behaviour was recorded 
for a wide range of (1q which makes our simulations reli- 
able in spite of the sign problem. In Fig.l (a-c) we present 
typical distributions of the first coefficients correspond- 
ing to n = 1 — 7 in cqn. (p"4|) (normalized with respect to 
the zeroth power coefficient) for (3 — 5 — 7.5 for differ- 
ent fxo- The coefficients are obtained using the bootstrap 
method on over 10000 configurations for /? = 5 increasing 
to over 30000 for (3 — 7.5. In spite of different values of 
the average sign in these simulations, the coefficients of 
the expansion (|l^ ) indicate good correspondence between 
coefficients obtained with different values of the update 
chemical potential /io: the normalized coefficients taken 
from different /iq values and equal power of the expansion 
variable correspond within the statistical error estimated 
using the bootstrap analysis. (To compare these coeffi- 
cients we had to shift the expansion by 2cosh/io/3.) 

We also performed a bootstrap analysis of the zeros in 
the /i plane which shows clear Gaussian distribution of 
their real and imaginary parts (see Fig. 2). In addition, 
we observe overlapping results (i.e. same zeros) obtained 
with different values of fiQ. The distribution of Yang- 
Lee zeros in the complex /i-plane is presented in Fig.3(a- 
c) for the zeros nearest to the real axis. We observe a 
gradual decrease of the imaginary part as the lattice size 
increases. The quantitative analysis of this behaviour is 
discussed in the next section. 



The critical domain can be identified by the behaviour 
of the density of Yang-Lee zeros' in the positive half-plane 
of the fugacity. We expect to find that this density is tem- 
perature and volume dependent as the system approaches 
the phase transition. If the temperature is much higher 
than the critical temperature, the zeros stay far from the 
positive real axis as it happens in the high-temperature 
limit of the one-dimensional Ising model (T c = 0) in 
which, for (3 — 0, the points of singularity of the free 
energy lie at fugacity value —1. As the temperature de- 
creases we expect the zeros to migrate to the positive 
half-plane with their density, in this region, increasing 
with the system's volume. 

Figures 4(a-c) show the number N{9) of zeros in the 
sector (0, 9) as a function of the angle 9. The zeros shown 
in these figures are those presented in Fig. 3 (a-c) in the 
chemical potential plane with other zeros lying further 
from the positive real half-axis added in . We included 
only the zeros having absolute value less than one which 
we are able to do because if yi is a zero in the fugacity 
plane, so is l/yi- The errors are shown where they were 
estimated using the bootstrap analysis (see Fig. 2). 

For (3 = 5, even for the largest simulated lattice 10 2 , 
all the zeros are in the negative half-plane. We notice a 
gradual movement of the pattern of the zeros towards the 
smaller 9 values with an increasing density of the zeros 
near 9 = ^ which is due to the increase of the lattice size. 
However, since no zeros are detected beyond this value 
of 9 we can assume that (3 — 5 is too high a temperature 
for a phase transition. 

At [3 — 6 we observe a qualitatively distinct behaviour 
for the two biggest lattice sizes (8 2 and 10 2 ). The low- 
lying zeros are in the positive half-plane with angular 
density estimated as the slope of the linear interpolation 
between the points on the graphs of N(9) for a given 
lattice size. These slopes increase with lattice size. Yet 
even the 10 2 result extrapolated towards 9 = yields a 
vanishing density of fugacity zeros on the positive real 
half- axis. 

The (3 — 7.5 results show the density of zeros increasing 
with the simulation volume at small values of 9. For the 
8 2 lattice the extrapolation of N(9) gives a non- vanishing 
slope at zero angle. We interpret these results as sig- 
nalling the existence of a criticality for this temperature. 
To obtain the scaling of the density of zeros near the pos- 
itive half-axis we would have to extend the simulation to 
larger lattices (work which is currently in progress). On 
the basis of the existing results, we conclude that the 
system described in our simulations approaches critical- 
ity for the temperature corresponding to the lowest of 
our simulations, namely (3 = 7.5. We do not exclude 
the possibility that larger lattice simulations may show a 
phase transition for an even slightly higher temperature, 
as the (3 = 6 results indicate. 
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IV. FINITE SIZE SCALING AND THE 
SINGULARITY OF THE DENSITY OF STATES 

As a starting point for the finite size analysis of the 
Yang-Lee singularities we recall the scaling hypothesis for 
the partition function singularities in the critical domain 
(□J. Following this hypothesis, for a change of scale of 
the linear dimension L 



we obtain a transformation of the partition function 

Zl — > Z L /\ 

which, in the critical domain, satisfies 

Z L {0,ii) = Z L/x (9 x ,fi x ) (17) 

where 9 = (ijr — 1), fl = (1 — jj-) and 9\ and jx x are 
the values of 9 and ft under the rescaling transforma- 
tion. This equation is correct up to a multiplication by 
an exponential stemming from the regular part of the free 
energy. In the scaling regime the following relationships 
hold ' 

6 X = \ ae 9 

A a = A a */J (18) 

A d -F s ing(0, A) = F s \ ng (9\,jl\) 

where ag and are the scaling exponents of the reduced 
temperature and reduced chemical potential respectively. 
While the A dependence of 9\ and jl\ expressed in the 
first two equations of the set (Q) are correct for A w 1, 
the third equation (in which -F s i ng is the singular part 
of the free energy density), follows in general from the 
conservation of the partition function under the scaling 
transformation. Here and in what follows we assume that 
this transformation does not generate new couplings in 
the model. Hence, the zeros of Zl will be the same as 
the zeros of Zl/ X - 

Assume that the simulation temperature is sufficiently 
low to allow for a phase transition with respect to the 
chemical potential. In this case, the zeros of the grand 
canonical partition function will correspond to the zeros 
of a scaling function Z 

Z L /x(9X a °,^') = Z(9L a \fiL a ») = 0. (19) 

Solution of this equation determines the relationship be- 
tween the reduced chemical potential and the tempera- 
ture through a function /$ 

jHL a » = h{6L a <>) (20) 

where i is the index of the root of the partition function. 
jjL c may vary with the temperature, and in particular we 



can regard /x c (T s i mu i a tion) as the critical chemical poten- 
tial at the temperature of simulation. Thus, for 9 = 
the previous equation becomes: 

JmL^ = MO) (21) 

which leads immediately to a determination of the expo- 
nent a^ from the logarithmic scaling of the zeros in the 
complex [i plane as 

\og\nL ~Mc|i = -"M lo g L + lo g/*(°)- ( 22 ) 

The notation of the last equation emphasizes the in- 
dependence of the chemical potential fi. 

The above derivation of the finite size scaling for the 
Yang-Lee zeros is a generalization of the ideas developed 
for a two-dimensional Ising model. In the present con- 
text, the essential difference between the Hubbard and 
Ising models, is in the phase boundary T c (/x). We do not 
expect this difference to significantly alter the scaling re- 
lations with respect to the reduced chemical potential. 
For the Ising model the phase boundary in the T — fi 
plane is a straight line fi — fi c — (with the identi- 
fication of fj, as the external magnetic field). This line 
extends from T = to T = T c , with the latter being the 
critical temperature of a second-order phase transition. 
Below T c the phase transition with respect to the exter- 
nal field is first-order accompanied by the characteristic 
two-phase coexistence. We expect that for the Hubbard 
model the boundary will be given by a line parametrised 
as T c (fi). As this line is not necessarily parallel to the 
T-axis, as one moves along either the T or fi directions, 
the phase boundary will be crossed. Since the free energy 
is an analytic function on both sides of the boundary, the 
limits of its derivatives are independent of the path along 
which these limits are taken (as long as the boundary is 
not crossed). Therefore, we would expect the same order 
of singularity and phase transition with respect to cither 
temperature or chemical potential. 

In what follows we assume a transition weaker than 
a first-order one (an assumption supported by simula- 
tion results of Refs.[H, |§ for the 2D Hubbard model) 
although even for a first-order phase transition the finite- 
size scaling of the zeros is expected to be correct. 

The critical exponent a M is related to the average par- 
ticle number via jl in the critical domain. We have to 
determine the relationship between a M and <5, where the 
latter is defined by 

< n >~ pi* (23) 

Following the real-space renormalization group treatment 
of Ref. |llj and assuming that the change of scale A is 
a continuous parameter, the exponent ag is related to 
the critical exponent v of the correlation length as ag = 
i . This result follows from the scaling of the correlation 
length £ as 
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(24) 



Thus, iterating the scaling of 9 (Eq. |l8|) and reaching the 
point A « 9~"o we obtain 



\6\-\ 



(25) 



(The absolute value in the last equation corresponds to 
the points above and below T c .) Now recalling Eq.jl^) 
and using the hyperscaling hypothesis (-F S ing ~ 9 vd ) we 
find: 



F^(9,fi) = \9\» d F sing (±l, 



\9\ va » 



(26) 



where 9\ has been scaled to ±1 and fl\ expressed in terms 
of // and 9. Differentiating this equation with respect to 
jl yields: 



< n > s 



|X=±1 



which determines the critical exponent of the parti- 
cle density with respect to the reduced temperature as 



v(d— a^). Substituting \9\ ~< n >^„g" into the ar- 
gument Y = g | ^ Q>J and demanding the hyperscaling 9- 
dependence to be preserved in cqn.(|26|) we obtain: 



< n >^ [i 



(28) 



which defines the critical exponent \ = —£- !L in terms of 
the scaling exponent of the Yang-Lee zeros. 

Fig. 5 presents the scaling of the imaginary part of the 
fj, zeros for different values of the temperature. The linear 
regression slope of the logarithm of the imaginary part of 
the zeros plotted against the logarithm of the inverse lin- 
ear dimension of the simulation volume, increases when 
the temperature decreases from fj = 5 to /3 = 6. The re- 
sults of (3 — 7.5 correspond to a M — 1 .3 within the errors 
of the zeros as the simulation volume increases from 6 2 
to 8 2 . As it is seen from Fig. 3, we can trace zeros with 
similar real part (Re(fii) « 0.7 which is also consistent 
with the critical value of the chemical potential given in 
Ref. |2^]) as the lattice size increases, which allows us to 
examine only the scaling of the imaginary part. Table 1 
presents the values of and i obtained from this figure. 



/3 




i 

/> 


5 


0.5 ±0.05 


3.0 ±0.4 


6 


1.3 ±0.2 


0.5 ±0.2 


7.5 


1.3 ±0.3 


0.5 ±0.4 



Table 1 The finite size scaling exponents and the critical 
exponent of particle-number singularity near criticality 
versus inverse temperature. 



We also note that the location of the zeros corresponds 
to the singularity of the electronic susceptibility \- I n 
Fig. 6 we show x, given by \ — d< ^]j > , as a function of 
the chemical potential on an 8 2 lattice. The location of 
the peaks of the susceptibility, rounded by the finite size 
effects, is in good agreement with the distribution of the 
real part of the Yang-Lee zeros in the complex //-plane 
(see Fig. 3) which is particularly evident in the (3 = 7.5 
simulations (Fig. 4(c)). The contribution of each zero to 
the susceptibility can be singled out by expressing the 
free energy as: 
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F= " Vi) 



(29) 



where y is the fugacity variable and yi is the correspond- 
ing zero of the partition function. The dotted lines on 
these plots correspond to the contribution of the nearby 
zeros while the full polynomial contribution is given by 
the solid lines. We see that the developing singularities 
are indeed governed by the zeros closest to the real axis. 
The sharpening of the singularity as the temperature de- 
creases is also in accordance with the dependence of the 
distribution of the zeros on the temperature. 

The singularities of the free energy and its derivative 
with respect to the chemical potential, can be related to 
the quasiparticle density of states. To do this we assume 
that single particle excitations accurately represent the 
spectrum of the system. The relationship between the 
average particle density and the density of states p{u>) is 
given by 



' " = / dui- ^7 tp(w) 



(30) 



which in the low-temperature regime [3 — > oo becomes 



< n >=< n > r 



eg 



< n > sing = 



du>p(u 



(31) 



Assuming < n > S i ng dependence upon the chemical po- 
tential as < 7i > s ing= (ft — Pc) 1 ^ one obtains the zero- 
temperature approximation for the density of states as 
the particle susceptibility x> namely 



Xsir 



1 



dfj, 



Psing(/j) OC -(n - H C )1 



(32) 



and hence the rate of divergence of the density of states. 

As in the case of Lifshitz transitions the singularity 
of the particle number is rounded at finite temperature. 
However, for sufficiently low temperatures, the singular- 
ity of the density of states remains manifest in the free 
energy, the average particle density, and particle suscep- 
tibility . The regular part of the density of states does 
not contribute to the criticality, so we can concentrate on 
the singular part only. Consider a behaviour of the type 
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PsingM oc (oj - fi c )* 1 for to > fi c and p s in g M = 
for uj < (j, c . The integral ( J30| ) will give in the region 
T << (fi — fi c ) the low-temperature leading order be- 
haviour for the particle density: 



< n > 



sing 



(33) 



with the value (5 for the particle number governed by the 
divergence of the density of states (at low temperatures) 
in spite of the finite-temperature rounding of the singu- 
larity itself. This rounding of the singularity is indeed 
reflected in the difference between the values of a M at 
(3 = 5 and = 6. 



V. DISCUSSION AND OUTLOOK 

We note that in our finite size scaling analysis we do 
not include logarithmic corrections. In particular, these 
corrections may prove significant when taking into ac- 
count the fact that we are dealing with a two-dimensional 
system in which the pattern of the phase transition is 
likely to be of Kosterlitz-Thouless type [^3). The loga- 
rithmic corrections to the scaling laws have been proven 
essential in a recent work of Kenna and Irving [^4| . In- 
clusion of these corrections would allow us to obtain the 
critical exponents with higher accuracy. However, such 
analysis would require simulations on even larger lattices. 

The linear fits for the logarithmic scaling and the criti- 
cal exponents obtained, are to be viewed as approximate 
values reflecting the general behaviour of the Yang-Lee 
zeros as the temperature and lattice size are varied. Al- 
though the bootstrap analysis provided us with accurate 
estimates of the statistical error on the values of the ex- 
pansion coefficients and the Yang-Lee zeros, the small 
number of zeros obtained with sufficient accuracy does 
not allow us to claim higher precision for the critical ex- 
ponents on the basis of more elaborate fittings of the scal- 
ing behaviour. The finite-size effects may still be signifi- 
cant, especially as the simulation temperature decreases, 
thus affecting the scaling of the Yang-Lee zeros with the 
system size. Larger lattice simulations will therefore be 
required for an accurate evaluation of the critical expo- 
nent for the particle density and the density of states. 
Nevertheless, the onset of a singularity at finite temper- 
ature, and its persistence as the lattice size increases, are 
evident. 

The estimate of the critical exponent for the diver- 
gence rate of the density of states of the quasiparticle 
excitation spectrum is particularly relevant to the high 
T c superconductivity scenario based on the van Hove sin- 
l7| . It is emphasized in Ref. 1 25 1 that 



gularities |2J,[f26| 
the logarithmic singularity of a two-dimensional electron 
gas can, due to electronic correlations, turn into a power- 
law divergence resulting in an extended saddle point at 
the lattice momenta (n,0) and (0,7r). In the case of the 



density of states diverging as the — ^ power, this singu- 
larity leads to a weak-BCS mechanism for superconduc- 
tivity with T c set at w 100K f||. The extended saddle 
point behaviour has been confirmed by experimental re- 
sults and numerical calculations 25 1, |^],Q,[^lJ. In 
numerical simulations the extended saddle point can be 
traced by a reconstruction of the quasiparticle spectral 
weight function from the simulation data on the Matsub- 
ara Green's functions. This procedure requires solving an 
ill-posed inverse problem (3^],|33| and is highly sensitive 
to the finite-doping simulation difficulties. The results 
obtained so far, cannot provide a quantitative descrip- 
tion of the singularity. The Yang-Lee analysis, however, 
provides a scheme for a direct evaluation of the critical 
exponents. 

We acknowledge stimulating discussions with A. 
Bratkovsky, R. Kenna, P. Kornilovitch and C. Creffield 
in the course of the work on this paper. This work was 
supported by SERC grant GR/J 18675. 
References 

1. R. Blankenbecler, D.J. Scalapino and R.L. Sugar, 
Phys.Rev. D40, 2278 (1981). 

2. J.E. Hirsch, Phys.Rev. B31, 4403 (1985). 

3. I.M. Barbour and E.G. Klepfish, Phys.Rev. B46, 
469 (1992). 

4. C.N. Yang and T.D. Lee, Phys. Rev. 87, 404 
(1952); 

T.D. Lee and C.N. Yang, Phys. Rev. 87, 410 
(1952). 

5. R.B. Griffiths, Jour. Math. Phys. 10, 1559 (1969). 

6. M. Suzuki and M.E. Fisher, Jour. Math. Phys. 12, 
235 (1971). 

7. F. Dunlop and CM. Newman, Commun. Math 
Phys. 44, 223 (1975). 

8. B. Simon and R.B. Griffiths, Commun. Math Phys. 
33, 145 (1973). 

9. K-C. Lee, Phys. Rev. Lett 73, 2801 (1994). 

10. N.D. Mcrmin and H. Wagner, Phys. Rev. Lett 17, 
1133 (1966). 

11. C. Itzykson, R.B. Pearson and J.B. Zuber, Nucl. 
Phys.B220, 415 (1983); 

C. Itzykson and J-M. Drouffe, Statistical Field The- 
ory v.l, Cambridge University Press 1991. 

12. M.L. Glasser, V. Privman and L.S. Schulman, 
Phys. Rev.B35 (1987) 1841 . 

13. R. Kenna and C.B. Lang, Nucl. Phys. B393, 461 
(1993). 



7 



14. I.M. Barbour, A.J. Bell and E.G. Klepfish, Nucl. 
Phys. B389, 285 (1993). 

15. I.M. Lifshitz, JETP 38, 1569 (1960). 

16. A. A. Abrikosov, Fundamentals of the Theory of 
Metals North-Holland (1988). 

17. P. Hall, The Bootstrap and Edgeworth expansion, 
Springer (1992). 

18. S.R. White et al, Phys.Rev. B40, 506 (1989). 

19. J.E. Hirsch, Phys. Rev. B28, 4059 (1983). 

20. M. Suzuki, Prog. Theor. Phys. 56, 1454 (1976). 

21. A. Moreo, D. Scalapino and E. Dagotto, Phys. 
i?eu.B43, 11442 (1991). 

22. N. Furukawa and M. Imada, J. Phys. Soc. Japan 
61, 3331 (1992). 

23. J. Kosterlitz and D. Thoulcss, J. Phys. C6, 1181 
(1973); 

J. Kosterlitz, J. Phys. C7, 1046 (1974). 

24. R. Kenna and A.C. Irving, unpublished. 

25. K. Gofron et al, Phys. Rev. Lett. 73, 3302 (1994). 

26. D.M. Newns, P.C. Pattnaik and C.C. Tsuei, Phys. 
Rev. B43, 3075 (1991); 

D.M. Newns et al, Phys. Rev. Lett. 24, 1264 

(1992) ; 

D.M. Newns et al, Phys. Rev. Lett. 73, 1264 
(1994). 

27. E. Dagotto, A. Nazarenko and A. Moreo, Phys. 
Rev. Lett. 74, 310 (1995). 

28. A. A. Abrikosov, J.C. Campuzano and K. Gofron, 
Physica (Amsterdam) 214C, 73 (1993). 

29. D.S. Dessau et al, Phys. Rev. Lett. 71, 2781 

(1993) ; 

D.M. King et al, Phys. Rev. Lett. 73, 3298 (1994); 
P. Aebi et al, Phys. Rev. Lett. 72, 2757 (1994). 

30. E. Dagotto, A. Nazarenko and M. Boninsegni, 
Phys. Rev. Lett. 73, 728 (1994). 

31. N. Bulut, D.J. Scalapino and S.R. White, Phys. 
Rev. Lett. 73, 748 (1994). 

32. S.R. White, Phys. Rev. B44, 4670 (1991); 

M. Vekic and S.R. White, Phys. Rev. B47, 1160 
(1993). 

33. C.E. Creffield, E.G. Klepfish, E.R. Pike and Sarben 
Sarkar, unpublished. 



Figure Captions 
Figure 1 

Bootstrap distribution of normalized coefficients for ex- 
pansion (|l4|) at different update chemical potential no for 
an 8 2 lattice. The corresponding power of expansion is 
indicated in the top figure, (a) f3 = 5, (b) (3 = 6, (c) 
13 = 7.5. 

Figure 2 

Bootstrap distributions for the Yang-Lee zeros in the 
complex n plane closest to the real axis, (a) 10 2 lat- 
tice at (3 = 5, (b) 10 2 lattice at (3 = 6, (c) 8 2 lattice at 
[3 = 7.5. 

Figure 3 

Yang-Lee zeros in the complex fj, plane closest to the real 
axis, (a) ,9 = 5, (b) (3 = 6, (c) (3 — 7.5. The correspond- 
ing lattice size is shown in the top right-hand corner. 

Figure 4 

Angular distribution of the Yang-Lee zeros in the com- 
plex fugacity plane Error bars are drawn where esti- 
mated, (a) (3 = 5, (b) [3 = 6, (c) [3 = 7.5. 

Figure 5 

Scaling of the imaginary part of fii (i?e(/ii) ks= 0.7) as 
a function of lattice size. a m u indicates the the fit of the 
logarithmic scaling. 

Figure 6 

Electronic susceptibility as a function of chemical poten- 
tial for an 8 2 lattice. The solid line represents the con- 
tribution of all the 2n x n y zeros and the dotted line the 
contribution of the six zeros nearest to the real-^ axis, 
(a) (3 = 5, (h) /3 = 6, (c) (3 = 7.5. 
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